# pragma once

# include <string>
# include <vector>
# include <cassert>

namespace linalg {

  class Array;
  class Matrix;

  namespace details {

    class IndexNumerator {
    public:

      IndexNumerator(size_t const nr, size_t const nc) throw()
        : rows_(nr), cols_(nc) {}

      size_t index(size_t const nr, size_t const nc) {
        return 0;
      }

    private:

      size_t rows_;
      size_t cols_;
    };



    /*! @brief Type of the function for index computation
     */
    typedef size_t (*indexer)(size_t, size_t, size_t, size_t);

    /*! @brief Computes index location in Fortran-like order
     */
    inline size_t fortranOrderIndex(size_t r, size_t c, size_t mr) throw () {
      assert(r < mr);

      return c * mr + r;
    }

    /*! @brief Computes item location in C-like order
     */
    inline size_t cOrderIndex(size_t r, size_t c, size_t mc) throw () {
      assert(c < mc);

      return r * mc + c;
    }

  }



  /*! @brief Array
   */
  class Array {
    friend class Matrix;
  public:

    enum class Mode {
      RAW,
      COLUMN
    };

  private:

    /*! @brief Array constructor
     */
    Array(double* ptr, Mode mode, size_t size, size_t Q) throw()
        : data_(ptr), size_(size), q_(Q), mode_(mode) {}

  public:

    /*! @brief Getter 
     */
    double operator()(size_t n) const throw() {
      if (mode_ == Mode::COLUMN)
        return data_[details::fortranOrderIndex(n, q_, size_)];
      else
        return data_[details::fortranOrderIndex(q_, n, size_)];
    }

    /*! @brief Setter
     */
    double& operator()(size_t n) throw() {
      if (mode_ == Mode::COLUMN)
        return data_[details::fortranOrderIndex(n, q_, size_)];
      else
        return data_[details::fortranOrderIndex(q_, n, size_)];
    }

  private:

    /*! @brief Pointer to data
     */
    double* const data_;

    /*! @brief Number of item (column or row)
     */
    size_t const size_;

    /*! @brief 
     */
    size_t const q_;
    

    /*! @brief Mode of array
     */
    Mode const mode_;
    
  };

  /*! @brief Matrix presentation
   */
  class Matrix {
    friend void svd(Matrix const&, Matrix&, Matrix&, Matrix&);
    friend void solve(Matrix const&, Matrix&);
  public:

    /*! @brief Default constructor
     */
    Matrix(): data_(), permutation_(), cols_(0), rows_(0) {}

    /*! @brief Matrix constructor
     *
     * @param[in] cols number of columns in the matrix
     * @param[in] rows number of rows in the matrix
     */
    Matrix(size_t cols, size_t rows) 
        : data_(cols*rows), permutation_(), cols_(cols), rows_(rows) {}

    /*! @brief Get the item of the matrix
     *
     * @param[in] i number of row 
     * @param[in] j number of column
     * @return reference to item of the matrix
     */
    inline double& operator()(size_t i, size_t j) throw() { 
      return data_[details::fortranOrderIndex(i, j, rows_)];
    }

    /*! @brief Get the item of the matrix
     *
     * @param[in] i number of row 
     * @param[in] j number of column
     * @return reference to item of the matrix
     */
    inline double operator()(size_t i, size_t j) const throw() { 
      return data_[details::fortranOrderIndex(i, j, rows_)];
    }

    /*! @brief Resize matrix
     *
     * @param[in] nr number of rows 
     * @param[in] nc number of columns
     */
    void resize(size_t const nr, size_t const nc);

    double const* data() const throw() { return &data_[0]; }
    double* data() throw() { return &data_[0]; }

    /*! @brief Get number of rows
     */
    inline size_t cols() const throw() { return cols_; }

    /*! @brief Get number of columns
     */
    inline size_t rows() const throw() { return rows_; }

    /*! @brief Compute LUD of the matrix
     */
    void LUD();

    /*! @brief Get column of the matrix
     *
     * @param[in] col number of column
     * @return column representation
     */
    Array column(size_t const col) throw () {
      return Array(&data_[0], Array::Mode::COLUMN, rows_, col);
    }

    /*! @brief Get row of the matrix
     *
     * @param[in] row number of row
     * @return row representation
     */
    Array row(size_t const row) throw() {
      return Array(&data_[0], Array::Mode::RAW, cols_, row);
    }

  private:

    //inline size_t index(size_t const r, size_t const c) const {
    //  assert(r < rows_);
    //  assert(c < cols_);
    //  return c*rows_+r; // Fortran mode
    //  return r*cols_+c; // C-mode
    //}

  private:

    std::vector<double> data_;
    std::vector<int> permutation_;

    size_t cols_;
    size_t rows_;
  };


  /*! @brief Compute SVD of the matrix
   *
   * @param[in] M matrix 
   * @param[out] u
   * @param[out] s singular values of matrix
   * @param[out] vt transpose version 
   */
  void svd(Matrix const& M, Matrix& u, Matrix& s, Matrix& vt);


  /*! @brief Solve linear system of equations (in place)
   *
   * @param[in] matrix matrix of the system
   * @param[in, out] right-hand side of the system
   */
  void solve(Matrix const& matrix, Matrix& rhs);
}
